#This computes rich subs moments and also other useful graphs can be restored
#Also outputs the  CV of x variables.
#
# Get BLP own elasticity average====
set.seed(seedn)
BLP.output <- monte.func.rich.subs.mlog(repi=1,year.monte=year.gph,objective="richsubs",vdgp=3)
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)]
BLP.priceadjust <- BLP.output.m[, mean(p_m * (1-s_m))] #  v=1 version --> wvp_m = 0 --> average own elasticity will be same as v=3
#
# coefficients of variation of the x 
CV <- function(x) round(sd(x)/mean(x),3)
CV.dt <- BLP.output.m[year==1990,lapply(.SD,CV),.SDcols=c("hpwt_m","air_m","mpd_m","space_m")]
setnames(CV.dt,c("hpwt_m","air_m","mpd_m","space_m"),c("hpwt","air","mpd","space"))
fwrite(CV.dt,"Tables/CF_BLPdata/CVx.csv")
#
#Generate the BLP richsubs moments for one year (last of their data) and 100 reps
#
#With full blown BLP heterogeneity
set.seed(seedn)
monte.func.rich.subs.mlog.moments(repi=1,year.monte=year.gph,vdgp = 3)
monte.allv <- monte.rep.rich.subs.moments(parallel=T,year.monte=year.gph,vdgp = 3)
monte.BLPmoments.v3 <- monte.allv[,.(crosselasmaha=round(mean(crosselasmaha),digits = 2) ,
                                     meanmaha=round(mean(meanmaha),digits = 2) ,
                                     IQR=round(mean(IQR),digits = 2),
                                     exvar=round(mean(exvar),digits = 3),
                                     covprob=round(mean(covprob),digits = 3),
                                     mncv=round(mean(mncv),digits = 3))]
#
#With only heterogeneity on demand for characteristics
set.seed(seedn)
monte.func.rich.subs.mlog.moments(repi=1,year.monte=year.gph,vdgp = 2)
monte.allv <- monte.rep.rich.subs.moments(parallel=T,year.monte=year.gph,vdgp = 2)
monte.BLPmoments.v2 <- monte.allv[,.(crosselasmaha=round(mean(crosselasmaha),digits = 2) ,
                                     meanmaha=round(mean(meanmaha),digits = 2) ,
                                     IQR=round(mean(IQR),digits = 2),
                                     exvar=round(mean(exvar),digits = 3),
                                     covprob=round(mean(covprob),digits = 3),
                                     mncv=round(mean(mncv),digits = 3))]
#
monte.func.rich.subs.mlog.moments(repi=1,year.monte=year.gph,vdgp = 1)
monte.BLPmoments.v1 <- data.table(crosselasmaha=0,meanmaha=monte.BLPmoments.v2$meanmaha,IQR=monte.BLPmoments.v2$IQR,exvar=0,covprob=0,mncv=0)
monte.BLPmoments <- data.table(v=1:3,rbind(monte.BLPmoments.v1,monte.BLPmoments.v2,monte.BLPmoments.v3))
fwrite(monte.BLPmoments,"Tables/CF_BLPdata/BLPrichsubs_moments.csv")
monte.BLPmoments
saveRDS(monte.BLPmoments,"Tables/CF_BLPdata/BLPrichsubs_moments.rds")
#
# Graphing the cross elas / maha regressions====
#
# V3 Log-linear
set.seed(seedn)
DT <- monte.func.rich.subs.mlog(repi=1,year.monte=year.gph,objective="richsubs",vdgp=3)
DTv3 <- copy(DT)
#regression of log cross elas against maha
res_loglin <- feols(log(cross_elas)~ maha_dist | m+j,DTv3[m !=j])  
summary(res_loglin) #fits better, more appropriate from theory? what are units?
#avplot regression by hand using fixest
res_1 <- feols(log(cross_elas)~ 1 | m+j,DTv3[m !=j])  
res_2 <- feols(maha_dist~ 1 | m+j,DTv3[m !=j])  
DT_res <- data.table(res_cross=res_1$residuals,res_maha=res_2$residuals,cross_elas = DTv3[m !=j]$cross_elas)
lmres2 <- lm(res_cross~ res_maha,DT_res)
lmres1 <- lm(log(cross_elas)~ res_maha,DT_res)
#avplot
pdf("Graphs/CF_BLPdata/richsubs_BLPdata_avplot_v3.pdf",height=5,width=5)
par(mar=c(4,5,1,1))  #bottom, left, top and right
DT_res[,plot(res_maha,res_cross,col=adjustcolor("grey",0.5),pch=20,
             ylab="log cross price elasticity (cond.)",xlab="Mahalanobis distance in attributes (cond.)")]
abline(lmres2$coefficients,col="black")
DT_res[,text(max(res_maha),max(res_cross),labels=paste("Coefficient: ",formatC(lmres2$coefficients[2],digits=3,format="f")),col="black",pos=2)]
dev.off()
#
#
# Graphing cross-elas of the v3 setting vs CES approx ====
#
# V3 graph
DTv3g <- copy(DTv3)
#
par(pty="s")
meanownelas_v3 <- DTv3g[ m==j , mean(own_elas_m)] # smaller than  BLP.meanownelas by construction (challenge 2)
DTv3g[,cross_elas := 100*cross_elas]
DTv3g[,cross_elas.ces := -100 * meanownelas_v3 * s_m]
pdf("Graphs/CF_BLPdata/richsubs_BLPdata_crosselas_v3.pdf",height=5,width=5)
par(mar=c(4,5,1,1))  #bottom, left, top and right
DTv3g[ m!=j ,HeadR::scatter(cross_elas.ces,cross_elas,col="grey",log="xy",
                            xlab="CES (log scale, X100)",
                            ylab=TeX("BLP (log scale, X100)"),
                            xaxt="n",yaxt="n")]
lab1 <-  c(0.01,0.1,1)
lab2 <- c(lab1,10)
axis(side =1,at = lab1,labels=lab1)
axis(side =2,at = lab2,labels=lab2,las=2)
abline(0,1)
dev.off()
#
#
